#### Claassen & Devine, The Politics of Imperial Nostalgia
#### Code to run main analyses 

library(car)
library(dplyr)
library(tidyr)
library(ggplot2)
library(forcats)
library(psych)
library(psychTools)
library(texreg)
library(arm)
library(lavaan)
library(semTable)
library(paran)
library(haven)
library(viridis)
library(RColorBrewer)
library(questionr)
library(list)
library(glmnet)
library(vip)
library(party)
library(permimp)
library(pdp)
library(ranger)
library(foreach)
library(doParallel)
library(tth)
library(plm)
library(clubSandwich)
library(stringr)
library(cjbart)
library(cregg)
library(ggpubr)
library(weights)
library(systemfit)
library(gridExtra)
library(patchwork) 
library(lavaan)
library(fastDummies)  
library(sandwich)
library(lmtest)
library(cowplot)


# folders

WD = rstudioapi::getActiveDocumentContext()$path 
setwd(dirname(WD))
print( getwd() )


## Read data

# read data
datm = read.csv("rude_w1-3_merg_scaled_rep.csv")

# reshape the dataset to long format
datl = datm %>%
  pivot_longer(cols = -ID, names_to = c(".value", "round"), 
               names_pattern = "(.*)_(w[1-3])") %>%
  mutate(round = as.numeric(gsub("w", "", round)))

datl$GOR = as.factor(datl$GOR)

# convert to panel dataframe
pdat = pdata.frame(datl, index = c("ID", "round"))



#### Party support correlations----

# rounds 1 & 2

vars1 = c("imp_att_w2", "left_ecoval_w2", "auth_val_w2", "immig_supp_w2", 
          "hos_sexism_w2", "popul_att_w2", "trust_pols_w2", "grp_polstat_w2", 
          "nat_id_eng_w2", "EU_indep_w2", "uk_ethant_7_w1")
vars2 = c("sup_Lab_w2", "sup_Cons_w2", "sup_Reform_w2", "sup_LibDem_w2", "sup_Green_w2")
cormat_w2 = data.frame(matrix(NA, nrow = length(vars1), ncol = length(vars2)))
for (i in 1:length(vars1)) {
  for (j in 1:length(vars2)) {
    cormat_w2[i, j] = round(wtd.cor(
      datm[[vars1[i]]], datm[[vars2[j]]], weight = datm$wt_w2)[1], 2)
  }
}
colnames(cormat_w2) =  c("Lab", "Cons", "Reform", "LibDem", "Green")
rownames(cormat_w2) = c(
  "imperial nostalgia", "left economic values", "authoritarian values", "immigration support", 
  "hostile sexism", "populist attitudes", "trust in politicians", "external political efficacy", 
  "English national identity", "independence from EU", "life better 50 years ago")
cormat_w2
write.csv(cormat_w2, "party_evals_corr_w2.csv")


# round 3

vars1 = c("imp_att_w3", "left_ecoval_w3", "auth_val_w3", "immig_supp_w3", 
          "nat_id_eng_w3", "nat_chauv_w3", "nat_pride_w3")
vars2 = c("sup_Lab_w3", "sup_Cons_w3", "sup_Reform_w3", "sup_LibDem_w3", "sup_Green_w3")
cormat_w3 = data.frame(matrix(NA, nrow = length(vars1), ncol = length(vars2)))
for (i in 1:length(vars1)) {
  for (j in 1:length(vars2)) {
    cormat_w3[i, j] = round(wtd.cor(
      datm[[vars1[i]]], datm[[vars2[j]]], weight = datm$wt_w3)[1], 2)
  }  
}
colnames(cormat_w3) =  c("Lab", "Cons", "Reform", "LibDem", "Green")
rownames(cormat_w3) = c(
  "imperial nostalgia", "left economic values", "authoritarian values", 
  "immigration support", "English national identity", "national chauvinism",
  "national pride")
cormat_w3
write.csv(cormat_w3, "party_evals_corr_w3.csv")


## plot combined correlations as heatmap

cormat_comb = rbind(cormat_w3, cormat_w2)
cormat_comb = cormat_comb[!grepl("1$", rownames(cormat_comb)), ]
cormat_comb$MAV = round(apply(cormat_comb, 1, function(x) mean(abs(x))), 2)
cormat_comb = cormat_comb[order(cormat_comb$MAV, decreasing = TRUE), ]
colnames(cormat_comb) =  c("Labour", "Conserv.", "Reform", "Lib. Dem.", "Green", "Mean")
cormat_comb
write.csv(cormat_comb, "party_evals_corr_w2&3.csv")

x_labs = c("Independence from\nEU", "Imperial nostalgia", "Immigration support", 
           "Authoritarian\nvalues", "Hostile sexism", "Left economic\nvalues", 
           "National chauvinism", "National pride", "English national\nidentity", 
           "Life better 50\nyears ago", "External efficacy", "Trust in\npoliticians", 
           "Populist attitudes")
plot_cols = colorRampPalette(RColorBrewer::brewer.pal(11, "BrBG"))(100)
plot_cols = viridis::viridis(100, option = "mako")[100:1]

# plot
png("cor_party_evals_heatmap.png", width = 1440, height = 1080)
par(mar = c(2, 9, 2, 1), cex = 2.4, tcl = -0.2)
image(x = 1:ncol(cormat_comb), y = 1:nrow(cormat_comb), 
      z = t(abs(cormat_comb))[, nrow(cormat_comb):1], 
      col = plot_cols, axes = FALSE, xlab = "", ylab = "",
      zlim = c(0, 1))
axis(3, at = 1:ncol(cormat_comb), labels = colnames(cormat_comb), tcl = 0, las = 1,
     mgp = c(1, 0.3, 0))
axis(2, at = 1:nrow(cormat_comb), labels = rev(x_labs), las = 2, 
     mgp = c(1, 0.4, 0))
mtext(side = 1, line = 0.4, cex = 2.4, 
      "Correlations with likelihood of ever voting for party")
for (i in 1:nrow(cormat_comb)) {
  for (j in 1:ncol(cormat_comb)) {
    text(j, nrow(cormat_comb) - i + 1, sprintf("%.2f", cormat_comb[i, j]), 
         cex = 0.9, col = "grey10")
  }
}
box()
dev.off()



#### Party support regressions----


## Wave 2 seemingly unrelated regressions

dat_w2 = datm[!(is.na(datm$wt_w2)), ]

# Create dummy variables for demogs
factor_vars = c("GOR_w2", "age_4cat_w2", "ethnic3_w2", "edu_3cat_w2", 
                 "religion_w2", "socgrade3_w2")
for(var in factor_vars) {
  if(var %in% names(dat_w2)) {
    dat_w2[[var]] = as.factor(dat_w2[[var]])
  }
}
all_dummies = c()
for(var in factor_vars) {
  var_levels = levels(dat_w2[[var]])[-1]
  for(level in var_levels) {
    clean_level = gsub("-", "_", level)
    clean_level = gsub(" ", "_", clean_level)
    clean_level = gsub("\\+", "_plus", clean_level)
    dummy_name = paste0(var, "_", clean_level)
    dat_w2[[dummy_name]] = as.numeric(dat_w2[[var]] == level)
    all_dummies = c(all_dummies, dummy_name)
  }
}
print(all_dummies)

predictor_string = paste("imp_att_w2 + left_ecoval_w2 + auth_val_w2 + immig_supp_w2 + 
                          hos_sexism_w2 + popul_att_w2 + nat_id_eng_w2 + EU_indep_w2 + 
                          political_attention_w2 + female_w2 + ownhome_w2 +", 
                          paste(all_dummies, collapse = " + "))

sur_model_syntax = paste0(
"sup_Cons_w2 ~ ", predictor_string, "
sup_Lab_w2 ~ ", predictor_string, "
sup_LibDem_w2 ~ ", predictor_string, "
sup_Reform_w2 ~ ", predictor_string, "
sup_Green_w2 ~ ", predictor_string, "

# Allow residual covariances (SUR structure)
sup_Cons_w2 ~~ sup_Lab_w2
sup_Cons_w2 ~~ sup_LibDem_w2
sup_Cons_w2 ~~ sup_Reform_w2
sup_Cons_w2 ~~ sup_Green_w2
sup_Lab_w2 ~~ sup_LibDem_w2
sup_Lab_w2 ~~ sup_Reform_w2
sup_Lab_w2 ~~ sup_Green_w2
sup_LibDem_w2 ~~ sup_Reform_w2
sup_LibDem_w2 ~~ sup_Green_w2
sup_Reform_w2 ~~ sup_Green_w2"
)

# fit SUR
sur_fit = sem(sur_model_syntax, data = dat_w2, sampling.weights = "wt_w2", 
               estimator = "MLR", missing = "fiml")
summary(sur_fit, fit.measures = TRUE, standardized = FALSE, rsquare = TRUE)

# Create table
extract.lavaan.sur = function(model, equation_name, include.rsquared = TRUE, ...) {
  params = parameterEstimates(model)
  reg_params = params[params$op == "~" & params$lhs == equation_name, ]
  coef_names = reg_params$rhs
  coef_vals = reg_params$est
  se_vals = reg_params$se
  pvals = reg_params$pvalue
  rsq_values = inspect(model, "r2")
  rsq = if(equation_name %in% names(rsq_values)) rsq_values[equation_name] else NA
  n_obs = lavInspect(model, "nobs")
  
  tr = createTexreg(
    coef.names = coef_names,
    coef = coef_vals,
    se = se_vals,
    pvalues = pvals,
    gof.names = c("$R^2$", "$N$"),
    gof = c(rsq, n_obs),
    gof.decimal = c(TRUE, FALSE)
  )
  return(tr)
}

sur_models = list(
  extract.lavaan.sur(sur_fit, "sup_Cons_w2"),
  extract.lavaan.sur(sur_fit, "sup_Lab_w2"),
  extract.lavaan.sur(sur_fit, "sup_LibDem_w2"),
  extract.lavaan.sur(sur_fit, "sup_Reform_w2"),
  extract.lavaan.sur(sur_fit, "sup_Green_w2")
)

custom_coef_names = c(
  "imp_att_w2" = "Imperial nostalgia",
  "left_ecoval_w2" = "Left-economic values", 
  "auth_val_w2" = "Authoritarian values",
  "immig_supp_w2" = "Immigration support",
  "hos_sexism_w2" = "Hostile sexism",
  "popul_att_w2" = "Populist attitudes",
  "nat_id_eng_w2" = "English identity",
  "EU_indep_w2" = "EU independence"
)

texreg(sur_models,
       file = "partypref_sur_w2.tex",
       custom.model.names = c("Conservative", "Labour", "Lib Dem", "Reform", "Green"),
       custom.coef.names = custom_coef_names,
       omit.coef = "GOR_w2_|socgrade3_w2_|age_4cat_w2_|edu_3cat_w2_|ethnic3_w2_|religion_w2_|ownhome_w2|female_w2|political_attention_w2",
       digits = 2,
       single.row = FALSE,
       leading.zero = FALSE,
       stars = c(0.05),
       dcolumn = TRUE,
       booktabs = TRUE,
       caption = "SUR Regression Results: Party Support",
       label = "tab:sur_results")



## Wave 2 SURs with nostalgic deprivation control

predictor_string = paste("imp_att_w2 + left_ecoval_w2 + auth_val_w2 + immig_supp_w2 + 
                          hos_sexism_w2 + popul_att_w2 + nat_id_eng_w2 + EU_indep_w2 + 
                          uk_ethant_7_w1 +political_attention_w2 + female_w2 + 
                          ownhome_w2 +", paste(all_dummies, collapse = " + "))

sur_model_syntax = paste0(
"sup_Cons_w2 ~ ", predictor_string, "
sup_Lab_w2 ~ ", predictor_string, "
sup_LibDem_w2 ~ ", predictor_string, "
sup_Reform_w2 ~ ", predictor_string, "
sup_Green_w2 ~ ", predictor_string, "

# Allow residual covariances (SUR structure)
sup_Cons_w2 ~~ sup_Lab_w2
sup_Cons_w2 ~~ sup_LibDem_w2
sup_Cons_w2 ~~ sup_Reform_w2
sup_Cons_w2 ~~ sup_Green_w2
sup_Lab_w2 ~~ sup_LibDem_w2
sup_Lab_w2 ~~ sup_Reform_w2
sup_Lab_w2 ~~ sup_Green_w2
sup_LibDem_w2 ~~ sup_Reform_w2
sup_LibDem_w2 ~~ sup_Green_w2
sup_Reform_w2 ~~ sup_Green_w2"
)

# fit SUR
sur_fit_2 = sem(sur_model_syntax, data = dat_w2, sampling.weights = "wt_w2", 
                 estimator = "MLR", missing = "fiml")
summary(sur_fit_2, fit.measures = TRUE, standardized = FALSE, rsquare = TRUE)

# Create table
sur_models_2 = list(
  extract.lavaan.sur(sur_fit_2, "sup_Cons_w2"),
  extract.lavaan.sur(sur_fit_2, "sup_Lab_w2"),
  extract.lavaan.sur(sur_fit_2, "sup_LibDem_w2"),
  extract.lavaan.sur(sur_fit_2, "sup_Reform_w2"),
  extract.lavaan.sur(sur_fit_2, "sup_Green_w2")
)

custom_coef_names = c(
  "imp_att_w2" = "Imperial nostalgia",
  "left_ecoval_w2" = "Left-economic values", 
  "auth_val_w2" = "Authoritarian values",
  "immig_supp_w2" = "Immigration support",
  "hos_sexism_w2" = "Hostile sexism",
  "popul_att_w2" = "Populist attitudes",
  "nat_id_eng_w2" = "English identity",
  "EU_indep_w2" = "EU independence",
  "uk_ethant_7_w1" = "Life better 50 years ago"
)

texreg(sur_models_2,
       file = "partypref_sur_2_w2.tex",
       custom.model.names = c("Conservative", "Labour", "Lib Dem", "Reform", "Green"),
       custom.coef.names = custom_coef_names,
       omit.coef = "GOR_w2_|socgrade3_w2_|age_4cat_w2_|edu_3cat_w2_|ethnic3_w2_|religion_w2_|ownhome_w2|female_w2|political_attention_w2",
       digits = 2,
       single.row = FALSE,
       leading.zero = FALSE,
       stars = c(0.05),
       dcolumn = TRUE,
       booktabs = TRUE,
       caption = "SUR Regression Results 2: Party Support",
       label = "tab:sur_results_2")


## wave 2 OLS regressions with SNP

datm$GOR_w2 = as.factor(datm$GOR_w2)
gen_mod_w2 = sup_Cons_w2 ~ imp_att_w2 + left_ecoval_w2 + auth_val_w2 + immig_supp_w2 + hos_sexism_w2 + popul_att_w2 + nat_id_eng_w2 + EU_indep_w2 + political_attention_w2 + age_4cat_w2 + ethnic3_w2 + female_w2 + edu_3cat_w2 + religion_w2 + socgrade3_w2 + ownhome_w2 + GOR_w2

vot_lm_con = lm(update(gen_mod_w2, sup_Cons_w2 ~ .), data = datm, weights = wt_w2)
vot_lm_lab = lm(update(gen_mod_w2, sup_Lab_w2 ~ .), data = datm, weights = wt_w2)
vot_lm_libd = lm(update(gen_mod_w2, sup_LibDem_w2 ~ .), data = datm, weights = wt_w2)
vot_lm_ref = lm(update(gen_mod_w2, sup_Reform_w2 ~ .), data = datm, weights = wt_w2)
vot_lm_grn = lm(update(gen_mod_w2, sup_Green_w2 ~ .), data = datm, weights = wt_w2)
vot_lm_snp = lm(update(gen_mod_w2, sup_SNP_w2 ~ . - GOR_w2), data = datm, weights = wt_w2)

vote_mods = list(vot_lm_con, vot_lm_lab, vot_lm_libd, vot_lm_ref, vot_lm_grn, vot_lm_snp)
screenreg(vote_mods, omit.coef = "GOR|socgrade|age_4cat|edu_3cat|ethnic|religion|ownhome|female",
          digits = 2, single.row = FALSE, leading.zero = FALSE)
texreg(file = "partypref_lm_w2.tex", vote_mods, digits = 2, leading.zero = FALSE,
       omit.coef = "GOR|socgrade|age_4cat|edu_3cat|ethnic|religion|ownhome|female|attention",
       custom.model.names = c("Cons", "Lab", "LDem", "Reform", "Green", "SNP"),
       custom.coef.names = c("Intercept", "Imperial nostalgia", "Left-economic values",
                             "Authoritarian values", "Immigration support",
                             "Hostile sexism", "Populist attitudes", "English identity",
                             "EU independence"),
       stars = (0.05), single.row = FALSE, dcolumn = TRUE, booktabs = TRUE)
htmlreg(file = "partypref_lm_w2.html", vote_mods, digits = 2, leading.zero = FALSE,
        omit.coef = "GOR|socgrade|age_4cat|edu_3cat|ethnic|religion|ownhome|female|attention",
        custom.model.names = c("Cons", "Lab", "LDem", "Reform", "Green", "SNP"),
        single.row = FALSE, stars = (0.05), caption = "Party pref OLS regression, w2")




#### Vote choice random forests----

## Actual party choice

# data from all waves (N = 663)

vars = c("imp_att_w3", "left_ecoval_w3", "auth_val_w3", "immig_supp_w3", 
         "nat_pride_w3", "nat_id_eng_w3", "nat_id_brit_w3", "nat_chauv_w3", 
         "political_attention_w3", "age_w3", "ethnic3_w3", "female_w3", 
         "edu_3cat_w3", "religion_w3", "socgrade6_w3", 
         "popul_att_w2", "trust_pols_w2", "hos_sexism_w2", "EU_indep_w2",
         "grp_polstat_w2", "uk_ethant_7_w1", "ownhome_w3", "GOR_w3", 
         "vote_intent_w3", "wt_w3")
rf_datm = datm[, vars]
rf_datm = rf_datm[complete.cases(rf_datm), ]
rf_datm = rf_datm %>% mutate(across(
  c(ethnic3_w3, socgrade6_w3, edu_3cat_w3, ownhome_w3, female_w3, 
    religion_w3, GOR_w3, vote_intent_w3), as.factor))
rf_datm = rf_datm %>% mutate(across(
  c(political_attention_w3, nat_id_eng_w3, nat_id_brit_w3, age_w3, 
    uk_ethant_7_w1, trust_pols_w2, EU_indep_w2), 
  as.numeric))
summary(rf_datm)
sapply(rf_datm, class)

# ranger model - probability / Brier score
rf_pref_int_w3 = ranger(vote_intent_w3 ~ . -wt_w3, rf_datm, num.trees = 500, replace = TRUE, 
                        case.weights = "wt_w3", importance = "permutation", probability = TRUE)
print(rf_pref_int_w3) # Brier score = 0.47

# extract vip scores
vip_choice_w3 = data.frame(var = names(rf_pref_int_w3$variable.importance), 
                           imp = as.numeric(rf_pref_int_w3$variable.importance))
vip_labs = c("Imperial nostalgia", "Left-economic values", 
             "Authoritarian values", "Immigration support", "National pride", 
             "English identity", "British identity", "Chauvinistic nationalism", 
             "Political attention", "Age", "Ethnicity", "Gender", 
             "Education level", "Religion", "Occupational grade", 
             "Populist attitudes", "Trust in politicians", "Hostile sexism",
             "EU independence", "Political efficacy", "Nostalgic deprivation", 
             "Homeownership", "Region")
vip_choice_w3$lab = vip_labs
vip_choice_w3 = vip_choice_w3[order(vip_choice_w3$imp, decreasing = TRUE), ]
vip_choice_w3 

# plot vip
vip_col = brewer.pal(5, "BrBG")[5]
png("vote_int_allwaves_vip_brier.png", height = 720, width = 720)
par(mfrow = c(1, 1), mar = c(3, 9, 2, 1.5), tcl = -0.2, cex = 1.9, las = 1, 
    mgp = c(1.8, 0.9, 0))
barplot(height = vip_choice_w3$imp[15:1], names.arg = vip_choice_w3$lab[15:1], 
        axes = FALSE, horiz = TRUE, space = 0.4, col = "white", 
        cex.names = 1, mgp = c(1, 0.2, 0), xlim = c(0, 0.04), ylim = c(0.7, 20.3))
abline(v = seq(0, 0.04, by = 0.01), lty = 3, lwd = 1, col = rgb(0.5, 0.5, 0.5, 1))
axis(side = 1, labels = TRUE, cex.axis = 1, mgp = c(1, 0.4, 0))
barplot(height = vip_choice_w3$imp[15:1], names.arg = NULL, axes = FALSE, horiz = TRUE, 
        space = 0.4, col = vip_col, add = TRUE)
mtext(side = 1, line = 1.5, "Variable importance for model's Brier score", cex = 1.9)
mtext(side = 3, line = 0.5, "Variables from all waves, N = 663", cex = 2.3)
dev.off()

# data from wave 3 (N = 1664)

vars = c("imp_att_w3", "left_ecoval_w3", "auth_val_w3", "immig_supp_w3", 
         "nat_pride_w3", "nat_id_eng_w3", "nat_id_brit_w3", "nat_chauv_w3", 
         "political_attention_w3", "age_w3", "ethnic3_w3", "female_w3", 
         "edu_3cat_w3", "religion_w3", "socgrade6_w3", "ownhome_w3", "GOR_w3", 
         "vote_intent_w3", "wt_w3")
rf_datm = datm[, vars]
rf_datm = rf_datm[complete.cases(rf_datm), ]
rf_datm = rf_datm %>% mutate(across(
  c(ethnic3_w3, socgrade6_w3, edu_3cat_w3, ownhome_w3, female_w3, 
    religion_w3, GOR_w3, vote_intent_w3), as.factor))
rf_datm = rf_datm %>% mutate(across(
  c(political_attention_w3, nat_id_eng_w3, nat_id_brit_w3, age_w3), 
  as.numeric))
summary(rf_datm)
sapply(rf_datm, class)

# ranger model - probability / Brier score
rf_pref_int_w3 = ranger(vote_intent_w3 ~ . -wt_w3, rf_datm, num.trees = 500, replace = TRUE, 
                        case.weights = "wt_w3", importance = "permutation", probability = TRUE)
print(rf_pref_int_w3) # Brier score = 0.47

# extract vip scores
vip_choice_w3 = data.frame(var = names(rf_pref_int_w3$variable.importance), 
                           imp = as.numeric(rf_pref_int_w3$variable.importance))
vip_choice_w3
vip_labs = c("Imperial nostalgia", "Left-economic values", 
             "Authoritarian values", "Immigration support", "National pride", 
             "English identity", "British identity", "Chauvinistic nationalism", 
             "Political attention", "Age", "Ethnicity", "Gender", 
             "Education level", "Religion", "Occupational grade", 
             "Homeownership", "Region")
vip_choice_w3$lab = vip_labs
vip_choice_w3 = vip_choice_w3[order(vip_choice_w3$imp, decreasing = TRUE), ]

# plot vip
vip_col = brewer.pal(5, "BrBG")[5]
png("vote_int_w3_vip_brier.png", height = 720, width = 720)
par(mfrow = c(1, 1), mar = c(3, 9, 2, 1.5), tcl = -0.2, cex = 1.9, las = 1, 
    mgp = c(1.8, 0.9, 0))
barplot(height = vip_choice_w3$imp[15:1], names.arg = vip_choice_w3$lab[15:1], 
        axes = FALSE, horiz = TRUE, space = 0.4, col = "white", 
        cex.names = 1, mgp = c(1, 0.2, 0), xlim = c(0, 0.04), ylim = c(0.7, 20.3))
abline(v = seq(0, 0.04, by = 0.01), lty = 3, lwd = 1, col = rgb(0.5, 0.5, 0.5, 1))
axis(side = 1, labels = TRUE, cex.axis = 1, mgp = c(1, 0.4, 0))
barplot(height = vip_choice_w3$imp[15:1], names.arg = NULL, axes = FALSE, horiz = TRUE, 
        space = 0.4, col = vip_col, add = TRUE)
mtext(side = 1, line = 1.5, "Variable importance for model's Brier score", cex = 1.9)
mtext(side = 3, line = 0.5, "Variables from 3rd wave only, N = 1664", cex = 2.3)
dev.off()



### Nostalgia correlates----


## Demographics

dem_vars = c("age_4cat", "gender", "edu_3cat", "GOR")
datl$gender = factor(ifelse(datl$female == 1, "Female", "Male"))
levels(datl$GOR) = c("North East", "North West", "Yorkshire & Humber",
                     "East Midlands", "West Midlands", "East of England",
                     "London", "South East", "South West", "Wales", "Scotland")
datl$edu_3cat = factor(str_to_title(datl$edu_3cat), 
                        levels = c("Low", "Medium", "High"))
datl$ethnic3 = factor(str_to_title(datl$ethnic3), 
                      levels = c("Other", "Other White", "White British"))
datl$age_4cat = factor(datl$age_4cat, levels = c("18-34", "35-49", "50-64", "65+"))
lapply(datl[, dem_vars], table)


# Pairwise combinations of demographic variables
pairs = combn(dem_vars, 2, simplify = FALSE)

# Map variable names to meaningful labels
label_map = c(
  "age_4cat" = "Age Group",
  "gender" = "Gender",
  "edu_3cat" = "Education Level",
  "GOR" = "Region"
)

# Create a list of plots
plot_list = lapply(pairs, function(pair) {

  # Fit the linear model
  formula = as.formula(paste0("imp_att_add ~ ", pair[1], "*", pair[2]))
  model_pair = lm(formula, data = datl)
  
  newdata_pair = expand.grid(
    var1 = unique(datl[[pair[1]]]),
    var2 = unique(datl[[pair[2]]])
  )
  newdata_pair = newdata_pair[complete.cases(newdata_pair), ]
  colnames(newdata_pair) = pair
  
  newdata_pair$predicted_effects = predict(model_pair, newdata_pair, type = "response")
  
  # Create the heatmap
  min_effect = 2
  max_effect = 4
  brbg_colors = brewer.pal(11, "BrBG")
  ggplot(newdata_pair, aes_string(x = pair[1], y = pair[2], fill = "predicted_effects")) +
    geom_tile() +
    scale_fill_gradientn(colors = brbg_colors, limits = c(min_effect, max_effect)) +
    labs(x = label_map[pair[1]], y = label_map[pair[2]], fill = "Predicted Effect") +
    theme_minimal() +  
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          axis.title = element_text(size = 14),
          axis.text = element_text(size = 13),
          plot.title = element_text(size = 16),
          legend.title = element_blank(),  
          axis.text.x = element_text(margin = margin(t = 0, unit = "pt")),
          axis.text.y = element_text(margin = margin(r = 0, unit = "pt"))) +
    ggtitle(paste(label_map[pair[1]], "vs.", label_map[pair[2]]))
})

# Remove NULL plots and arrange the remaining plots
png("imp_att_by_2way_demogs.png", width = 1080, height = 1080, 
    pointsize = 16)
plot_list = plot_list[!sapply(plot_list, is.null)]
grid.arrange(grobs = plot_list, ncol = 2)
dev.off()


## Nostalgia by demographics dot plot, wave 2

dem_vars = c("age_4cat_w2", "GOR_w2", "ethnic3_w2", "gender_w2", "edu_3cat_w2")
datm$gender_w2 = factor(ifelse(datm$female_w2 == 1, "Female", "Male"))
datm$GOR_w2 = as.factor(datm$GOR_w2)
levels(datm$GOR_w2) = c("North East", "North West", "Yorks & Humb.",
                     "East Midlands", "West Midlands", "East England",
                     "London", "South East", "South West", "Wales", "Scotland")
datm$edu_3cat_w2 = factor(str_to_title(datm$edu_3cat_w2), 
                       levels = c("Low", "Medium", "High"))
datm$ethnic3_w2 = factor(str_to_title(datm$ethnic3_w2), 
                      levels = c("Other", "Other White", "White British"))
datm$age_4cat_w2 = factor(datm$age_4cat_w2, levels = c("18-34", "35-49", "50-64", "65+"))
lapply(datm[, dem_vars], table)

demog_mods = list()
for (var in dem_vars) {
  if (!is.factor(datm[[var]])) {
    datm[[var]] = as.factor(datm[[var]])
  }
  mod_base = lm(imp_att_add_w2 ~ -1 + get(var), data = datm, weights = wt_w2)
  demog_mods[[var]] = list(coefs = coef(mod_base), ses = sqrt(diag(vcov(mod_base))))
}
demog_mods

# plot

png("imp_att_add_by_demogs_w2.png", width = 1080, height = 720)
plot_col = brewer.pal(5, "BrBG")[5]
trans_plot_col = adjustcolor(plot_col, alpha.f = 0.3)  
titles = c("Age", "Region", "Ethnicity", "Gender", "Education")
layout(matrix(c(1, 2, 2, 3, 4, 5), ncol = 2, byrow = FALSE))
for (i in 1:length(dem_vars)) {
  par(mar = c(1.5, 5, 1.5, .5), tcl = -0.2, cex = 1.8, las = 1)
  coefs = demog_mods[[dem_vars[i]]]$coefs
  ses = demog_mods[[dem_vars[i]]]$ses
  inds = 1:length(demog_mods[[i]][[1]])
  data_viol = datm[complete.cases(datm[, c("imp_att_add_w2", dem_vars[i])]), ]
  plot(coefs, inds, pch = 16, axes = FALSE, col = "white",
       xlab = "", ylab = "", main = NULL, xaxs = "i",
       xlim = c(1, 5), ylim = c(0.5, length(inds)+0.5), adj = 1)
  abline(v = 2:4, lty = 3, lwd = 1, col = rgb(.7, .7, .7, 1))
  if (nrow(data_viol) > 0) {
    for (j in 1:length(levels(data_viol[[dem_vars[i]]]))) {
      subgrp_data = data_viol[data_viol[[dem_vars[i]]] == 
                                  levels(data_viol[[dem_vars[i]]])[j], ]
      subgrp_wt = subgrp_data$wt_w2 / sum(subgrp_data$wt_w2, na.rm = TRUE)
      wt_dens = density(subgrp_data$imp_att_add_w2, 
                        weights = subgrp_wt, adjust = 0.5)
      scale_dens = wt_dens$y / max(wt_dens$y)
      polygon(c(wt_dens$x, rev(wt_dens$x)),
              c((j + (1 - scale_dens * 0.5) - 1), rev(j + scale_dens * 0.5)), 
              col = trans_plot_col, border = NA)
    }
  }
  points(coefs, inds, pch = 16, cex = 1.1, col = plot_col)
  segments(coefs - 1.96*ses, 1:length(coefs), coefs + 1.96*ses, 1:length(coefs), 
           lwd = 2, col = plot_col)
  mtext(titles[i], side = 3, line = 0.2, cex = 1.8, adj = 0.5)
  axis(side = 1, at = 1:5, labels = TRUE, cex.axis = 0.9, mgp = c(1.5, 0.2, 0))
  axis(side = 2, at = inds, labels = levels(datm[[dem_vars[i]]])[inds], 
       cex.axis = 0.9, lwd = 0, lwd.ticks = 0, mgp = c(1, 0.1, 0))
  box()
}
dev.off()




## Attitudinal variables

# round 2

covars_w2 = c("imp_att_w2", "imp_emot_diff_w2", "trust_pols_w2", "satis_dem_w2", 
              "dem_sup_w2", "auth_val_w2", "immig_supp_w2", "left_ecoval_w2", 
              "grp_polstat_w2", "hos_sexism_w2", "popul_att_w2", "nat_id_eng_w2", 
              "nat_id_brit_w2", "EU_indep_w2", "uk_ethant_7_w1")

nost_variables = covars_w2[1:2]
corlist_w2 = list()
for (n_var in nost_variables) {
  cor_vector = numeric(length = length(covars_w2) - 2) 
  for (j in 3:length(covars_w2)) {
    x = datm[[n_var]]
    y = datm[[covars_w2[j]]]
    w = datm$wt_w2
    cor_vector[j - 2] = wtd.cor(x, y, weight = w)[1]
  }
  corlist_w2[[n_var]] = cor_vector
}
cors_w2 = as.matrix(as.data.frame(corlist_w2))

# round 3

covars_w3 = c("imp_att_w3", "imp_emot_diff_w3", "nat_chauv_w3", "nat_pride_w3")

nost_variables = covars_w3[1:2]
corlist_w3 = list()
for (n_var in nost_variables) {
  cor_vector = numeric(length = length(covars_w3) - 2) 
  for (j in 3:length(covars_w3)) {
    x = datm[[n_var]]
    y = datm[[covars_w3[j]]]
    w = datm$wt_w3
    cor_vector[j - 2] = wtd.cor(x, y, weight = w)[1]
  }
  corlist_w3[[n_var]] = cor_vector
}
cors_w3 = as.matrix(as.data.frame(corlist_w3))

# combine, label and sort
cor_mat = rbind(cors_w2, cors_w3) %>%
  `rownames<-`(c(covars_w2[3:15], covars_w3[3:4])) %>%  
  `colnames<-`(c("imp_att", "imp_emot")) %>% 
  as.data.frame() %>%  
  arrange(desc(imp_att)) 



## plot correlations as heatmap

var_labs = c("Authoritarian values", "National chauvinism", "EU independence",
             "National pride", "Hostile sexism", "British identity", 
             "English identity", "Satisfaction with\ndemocracy", 
             "Life better 50 years\nago", "Trust in politicians", "External efficacy",
             "Populist attitudes", "Democratic support", "Left-economic values", 
             "Immigration support")
cor_mat_abs = cor_mat[order(abs(cor_mat$imp_att), decreasing = TRUE), ]
var_labs_abs = var_labs[order(abs(cor_mat$imp_att), decreasing = TRUE)]
plot_cols = viridis::viridis(100, option = "mako")[100:1]

png("nostalgia_cor_cov_heatmap.png", width = 1080, height = 1080)
par(mar = c(2, 11, 2.5, 1), cex = 2.4, tcl = -0.2)
image(x = 1:ncol(cor_mat_abs), y = 1:nrow(cor_mat_abs), 
      z = t(abs(cor_mat_abs))[, nrow(cor_mat_abs):1], 
      col = plot_cols, axes = FALSE, xlab = "", ylab = "",
      zlim = c(0, 1))
axis(3, at = 1:ncol(cor_mat_abs), labels = c("Imperial\nAttitudes", "Imperial\nEmotions"), 
     tcl = 0, las = 1, mgp = c(1.4, 0.3, 0))
axis(2, at = 1:nrow(cor_mat_abs), labels = rev(var_labs_abs), las = 2, 
     mgp = c(1, 0.4, 0), cex.axis = 0.9)
for (i in 1:nrow(cor_mat_abs)) {
  for (j in 1:ncol(cor_mat_abs)) {
    text(j, nrow(cor_mat_abs) - i + 1, sprintf("%.2f", cor_mat_abs[i, j]), 
         cex = 0.9, col = "grey10")
  }
}
box()
dev.off()



### Conjoint main effect analysis----


## Clean and recode


dat3 = datm %>% filter(part_w3 == 1)

# select conjoint data 
cnj_covars = c("imp_att_w3", "pref_party_w3", "vote_intent_w3",
               "sup_Cons_w3", "sup_Lab_w3", "sup_Reform_w3", 
               "sup_LibDem_w3", "sup_Green_w3", "auth_val_w3", 
               "nat_chauv_w3", "imp_emot_neg_w3", 
               "imp_emot_pos_w3", "wt_w3")
cnj_vars = c("q_attr1_concept1_task1_w3", "q_attr1_concept1_task2_w3",
             "q_attr1_concept1_task3_w3", "q_attr1_concept2_task1_w3", 
             "q_attr1_concept2_task2_w3", "q_attr1_concept2_task3_w3",
             "q_attr2_concept1_task1_w3", "q_attr2_concept1_task2_w3",
             "q_attr2_concept1_task3_w3", "q_attr2_concept2_task1_w3",
             "q_attr2_concept2_task2_w3", "q_attr2_concept2_task3_w3", 
             "q_attr3_concept1_task1_w3", "q_attr3_concept1_task2_w3",
             "q_attr3_concept1_task3_w3", "q_attr3_concept2_task1_w3", 
             "q_attr3_concept2_task2_w3", "q_attr3_concept2_task3_w3", 
             "q_attr4_concept1_task1_w3", "q_attr4_concept1_task2_w3",
             "q_attr4_concept1_task3_w3", "q_attr4_concept2_task1_w3",
             "q_attr4_concept2_task2_w3", "q_attr4_concept2_task3_w3", 
             "q_attr5_concept1_task1_w3", "q_attr5_concept1_task2_w3",
             "q_attr5_concept1_task3_w3", "q_attr5_concept2_task1_w3", 
             "q_attr5_concept2_task2_w3", "q_attr5_concept2_task3_w3", 
             "q_attr6_concept1_task1_w3", "q_attr6_concept1_task2_w3",
             "q_attr6_concept1_task3_w3", "q_attr6_concept2_task1_w3", 
             "q_attr6_concept2_task2_w3", "q_attr6_concept2_task3_w3", 
             "q_attr7_concept1_task1_w3", "q_attr7_concept1_task2_w3", 
             "q_attr7_concept1_task3_w3", "q_attr7_concept2_task1_w3", 
             "q_attr7_concept2_task2_w3", "q_attr7_concept2_task3_w3", 
             "q_attr8_concept1_task1_w3", "q_attr8_concept1_task2_w3",
             "q_attr8_concept1_task3_w3", "q_attr8_concept2_task1_w3", 
             "q_attr8_concept2_task2_w3", "q_attr8_concept2_task3_w3", 
             "Q1_taskeen1_w3",            "Q1_taskeen1_1_w3",
             "Q1_taskeen1_2_w3",          "Q2_taskeen1_w3",            
             "Q1_taskeen2_w3",            "Q1_taskeen2_1_w3",          
             "Q1_taskeen2_2_w3",          "Q2_taskeen2_w3",           
             "Q1_taskeen3_w3",            "Q1_taskeen3_1_w3",          
             "Q1_taskeen3_2_w3",        "Q2_taskeen3_w3")

dat_cnj1 = dat3[, c("ID", cnj_covars, cnj_vars)] 
names(dat_cnj1)

# rename the dvs
dv_nam1 = c("Q1_taskeen1_1_w3", "Q1_taskeen2_1_w3", "Q1_taskeen3_1_w3", "Q1_taskeen1_2_w3", 
            "Q1_taskeen2_2_w3", "Q1_taskeen3_2_w3", "Q2_taskeen1_w3", "Q2_taskeen2_w3", 
            "Q2_taskeen3_w3")
dv_nam2 = c("CExp_rating_MP1_task1", "CExp_rating_MP1_task2", "CExp_rating_MP1_task3", 
            "CExp_rating_MP2_task1", "CExp_rating_MP2_task2", "CExp_rating_MP2_task3", 
            "CExp_MPchoice_task1", "CExp_MPchoice_task2", "CExp_MPchoice_task3")
names_vector = setNames(dv_nam2, dv_nam1)
dat_cnj1 <- dat_cnj1 %>%
  dplyr::rename_with(~ names_vector[.x], .cols = dplyr::all_of(dv_nam1))
names(dat_cnj1)

# remove wave identifier from the vars
colnames(dat_cnj1) = gsub("_w3$", "", colnames(dat_cnj1))
names(dat_cnj1)

# pivot to long
cnj_long = dat_cnj1 %>%
  pivot_longer(
    cols = -c(ID, wt, imp_att, auth_val, nat_chauv, pref_party, vote_intent, 
              sup_Cons, sup_Lab, sup_Reform, sup_LibDem, sup_Green, 
              imp_emot_neg, imp_emot_pos),
    names_to = c(".value", "replicate"), 
    names_pattern = "(.*)(_task\\d)$" 
  ) %>%
  arrange(ID, replicate)

names(cnj_long) = gsub("concept", "MP", names(cnj_long))
names(cnj_long)

# separate the MP1 and MP2 attributes into two datasets
mp1 = cnj_long %>%
  dplyr::select(ID, wt, imp_att, auth_val, nat_chauv, pref_party, vote_intent, 
                sup_Cons, sup_Lab, sup_Reform, sup_LibDem, sup_Green, 
                imp_emot_neg, imp_emot_pos, CExp_MPchoice, 
                replicate, ends_with("_MP1")) %>%
  mutate(MP_profile = 1)

mp2 = cnj_long %>%
  dplyr::select(ID, wt, CExp_MPchoice, imp_att, auth_val, nat_chauv, pref_party, 
                vote_intent, sup_Cons, sup_Lab, sup_Reform, sup_LibDem, sup_Green, 
                imp_emot_neg, imp_emot_pos, replicate, 
                ends_with("_MP2")) %>%
  mutate(MP_profile = 2)

# rename the columns to remove the MP suffix
mp1 = rename_at(mp1, vars(ends_with("_MP1")), ~str_replace(., "_MP1", ""))
mp2 = rename_at(mp2, vars(ends_with("_MP2")), ~str_replace(., "_MP2", ""))

# combine the datasets
dat_cnj = bind_rows(mp1, mp2)

# recode the replicate variable
dat_cnj = dat_cnj %>%
  mutate(replicate = as.integer(str_replace(replicate, "_task", "")))

# create the MP_chosen variable
dat_cnj = dat_cnj %>%
  mutate(MP_chosen = as.integer(CExp_MPchoice == MP_profile))

# drop the CExp_MPchoice column
dat_cnj = dplyr::select(dat_cnj, -CExp_MPchoice)

# recode the attributes
dat_cnj = dat_cnj %>%
  mutate(
    MP_educ = dplyr::recode(q_attr1, `1` = "Left school 16", `2` = "College", 
                            `3` = "University"),
    MP_occup = dplyr::recode(q_attr2, `1` = "Doctor", `2` = "Accountant", `3` = "Chef", 
                             `4` = "Supermarket cashier"),
    MP_party = dplyr::recode(q_attr3, `1` = "Labour", `2` = "Conservative", 
                             `3` = "Liberal Democrat"),
    MP_ethnic = dplyr::recode(q_attr4, `1` = "White", `2` = "Black", `3` = "Asian"),
    MP_gender = dplyr::recode(q_attr5, `1` = "Female", `2` = "Male"),
    MP_age = dplyr::recode(q_attr6, `1` = "27", `2` = "43", `3` = "57", `4` = "72"),
    MP_tax_views = dplyr::recode(q_attr7, `1` = "Increase spending", 
                                 `2` = "Balance tax and spend", `3` = "Cut taxes"),
    MP_empire_views = dplyr::recode(q_attr8, `1` = "Atrocities", `2` = "Both", 
                                    `3` = "Civilising effect")
  )

# Drop the old attribute columns
dat_cnj = as.data.frame(dplyr::select(dat_cnj, -starts_with("q_attr")))

# Factorise and relevel attributes
attr_var = c("MP_educ", "MP_occup", "MP_party", "MP_ethnic", "MP_gender", 
             "MP_age", "MP_tax_views", "MP_empire_views")
dat_cnj = dat_cnj %>%
  mutate(across(all_of(attr_var), factor))
dat_cnj$MP_ethnic = relevel(dat_cnj$MP_ethnic, ref = "White")
dat_cnj$MP_empire_views = relevel(dat_cnj$MP_empire_views, ref = "Both")

# drop profiles with duplicate values
dat_cnj = dat_cnj %>% filter(!is.na(replicate))

# drop profiles with DKs for ratings
dat_cnj$CExp_rating[dat_cnj$CExp_rating == 99] = NA
dat_cnj_tr = dat_cnj %>% filter(!is.na(CExp_rating))

# Sort the datasets by ID, then by replicate, then by MP_profile
dat_cnj = dat_cnj %>%
  arrange(ID, replicate, MP_profile) %>%
  mutate(
    MP_empire_views = factor(MP_empire_views, 
                             levels = c("Civilising effect", "Both", "Atrocities")),
    MP_educ = factor(MP_educ, 
                     levels = c("University", "College", "Left school 16")),
    MP_occup = factor(MP_occup, 
                      levels = c("Doctor", "Accountant", "Chef", "Supermarket cashier")),
    MP_tax_views = factor(MP_tax_views, 
                          levels = c("Increase spending", "Balance tax and spend", "Cut taxes"))
  )

dat_cnj_tr = dat_cnj_tr %>%
  arrange(ID, replicate, MP_profile) %>%
  mutate(
    MP_empire_views = factor(MP_empire_views, 
                             levels = c("Civilising effect", "Both", "Atrocities")),
    MP_educ = factor(MP_educ, 
                     levels = c("University", "College", "Left school 16")),
    MP_occup = factor(MP_occup, 
                      levels = c("Doctor", "Accountant", "Chef", "Supermarket cashier")),
    MP_tax_views = factor(MP_tax_views, 
                          levels = c("Increase spending", "Balance tax and spend", "Cut taxes"))
  )




## Forced choice, marginal means, full dataset

plot_col = brewer.pal(5, "BrBG")[5]

mm1 = cj(MP_chosen ~ MP_empire_views + MP_educ + MP_occup + MP_party + MP_ethnic + MP_gender
          + MP_age + MP_tax_views, 
          data = dat_cnj, id = ~ID, estimate = "mm", weights = ~wt,
          feature_labels = list("MP_empire_views" = "Empire\nviews",
                                "MP_educ" = "Education", 
                                "MP_occup" = "Occupation",
                                "MP_party" = "Party",
                                "MP_ethnic" = "Ethnicity",
                                "MP_gender" = "Gender", 
                                "MP_age" = "Age",
                                "MP_tax_views" = "Tax\nviews"))

conjfig1.1 = mm1 %>% 
  ggplot(aes(level)) +
  geom_vline(xintercept = 0.5, color="grey10") +
  geom_segment(aes(x = lower, xend = upper, y = level, yend = level), 
               color = ifelse(mm1$feature %in% "Empire\nviews", plot_col, "grey30"), 
               size = 2.5, alpha = .8) +
  geom_point(aes(x = estimate, y = level), 
             color = ifelse(mm1$feature == "Empire\nviews", plot_col, "grey30"), 
             fill = "white", size = 2.5, pch = 21) +
  facet_grid(feature~., scales='free', space='free') +
  labs(y = NULL, x ='Marginal means', title = "Forced choice") +
  theme_classic() +
  theme(
    strip.placement = 'outside',
    panel.background = element_rect(color='black'),
    legend.position = "bottom",
    legend.title = element_blank(),
    strip.background= element_blank(),
    axis.text.x = element_text(size = 10),
    strip.text = element_blank(),
    axis.title.x = element_text(face="bold", size = 10),
    axis.text.y = element_text(face="bold", size = 10)) +
  theme(legend.position = "none")

ggsave("conj_choice.png", 
       dpi=500, height=20, width=20, units="cm")


## Ratings, marginal means, trimmed dataset

mm2 = cj(CExp_rating ~ MP_empire_views + MP_educ + MP_occup + MP_party + MP_ethnic + MP_gender
          + MP_age + MP_tax_views, 
          data = dat_cnj_tr, id = ~ID, estimate = "mm", 
          feature_labels = list("MP_empire_views" = "Empire\nviews",
                                "MP_educ" = "Education", 
                                "MP_occup" = "Occupation",
                                "MP_party" = "Party",
                                "MP_ethnic" = "Ethnicity",
                                "MP_gender" = "Gender", 
                                "MP_age" = "Age",
                                "MP_tax_views" = "Tax\nviews"))

conjfig1.2 = mm2 %>% 
  ggplot(aes(level)) +
  geom_vline(xintercept = 5.83, color="grey10") +
  geom_segment(aes(x = lower, xend = upper, y = level, yend = level), 
               color = ifelse(mm2$feature %in% "Empire\nviews", plot_col, "grey30"), 
               size = 2.5, alpha = .8) +
  geom_point(aes(x = estimate, y = level), 
             color = ifelse(mm2$feature == "Empire\nviews", plot_col, "grey30"), 
             fill = "white", size = 2.5, pch = 21) +
  facet_grid(feature ~ ., scales='free', space='free') +
  labs(y=NULL, x='Marginal means', title = "Profile rating") +
  theme_classic() +
  theme(
    strip.placement = 'outside',
    panel.background = element_rect(color='black'),
    legend.position = "bottom",
    legend.title = element_blank(),
    strip.background = element_blank(),
    strip.text = element_text(size = 12),
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(face="bold", size = 10),
    axis.ticks.y = element_blank(),
    axis.text.y = element_blank()) +
  theme(legend.position = "none")

ggsave("conj_rating.png", 
       dpi=500, height=20, width=25, units="cm")


## Combining the choice and the ranking

fullchoice = rbind(mm1, mm2)
ggarrange(conjfig1.1, 
          conjfig1.2,
          widths = c(1, 0.8))
ggsave("conj_comb.png", 
       dpi=500, height=20, width=20, units="cm")


## Forced choice, AMCEs, full dataset

plot_col = brewer.pal(5, "BrBG")[5]

amce1 = cj(MP_chosen ~ MP_empire_views + MP_educ + MP_occup + MP_party + MP_ethnic + MP_gender
            + MP_age + MP_tax_views, 
            data = dat_cnj, id = ~ID, estimate = "amce", weights = ~wt,
            feature_labels = list("MP_empire_views" = "Empire\nviews",
                                "MP_educ" = "Education", 
                                "MP_occup" = "Occupation",
                                "MP_party" = "Party",
                                "MP_ethnic" = "Ethnicity",
                                "MP_gender" = "Gender", 
                                "MP_age" = "Age",
                                "MP_tax_views" = "Tax\nviews"))

conjfig1.1 = amce1 %>% 
  ggplot(aes(level)) +
  geom_vline(xintercept = 0, color="grey10") +
  geom_segment(aes(x = lower, xend = upper, y = level, yend = level), 
               color = ifelse(amce1$feature %in% "Empire\nviews", plot_col, "grey30"), 
               size = 2.5, alpha = .8) +
  geom_point(aes(x = estimate, y = level), 
             color = ifelse(amce1$feature == "Empire\nviews", plot_col, "grey30"), 
             fill = "white", size = 2.5, pch = 21) +
  facet_grid(feature~., scales='free', space='free') +
  labs(y = NULL, x ='AMCE', title = "Forced choice") +
  theme_classic() +
  theme(
    strip.placement = 'outside',
    panel.background = element_rect(color='black'),
    legend.position = "bottom",
    legend.title = element_blank(),
    strip.background= element_blank(),
    axis.text.x = element_text(size = 10),
    strip.text = element_blank(),
    axis.title.x = element_text(face="bold", size = 10),
    axis.text.y = element_text(face="bold", size = 10)) +
  theme(legend.position = "none")


## Ratings, AMCEs, trimmed dataset

amce2 = cj(CExp_rating ~ MP_empire_views + MP_educ + MP_occup + MP_party + MP_ethnic + MP_gender
            + MP_age + MP_tax_views, 
            data = dat_cnj_tr, id = ~ID, estimate = "amce", 
            feature_labels = list("MP_empire_views" = "Empire\nviews",
                                "MP_educ" = "Education", 
                                "MP_occup" = "Occupation",
                                "MP_party" = "Party",
                                "MP_ethnic" = "Ethnicity",
                                "MP_gender" = "Gender", 
                                "MP_age" = "Age",
                                "MP_tax_views" = "Tax\nviews"))

conjfig1.2 = amce2 %>% 
  ggplot(aes(level)) +
  geom_vline(xintercept = 0, color="grey10") +
  geom_segment(aes(x = lower, xend = upper, y = level, yend = level), 
               color = ifelse(amce2$feature %in% "Empire\nviews", plot_col, "grey30"), 
               size = 2.5, alpha = .8) +
  geom_point(aes(x = estimate, y = level), 
             color = ifelse(amce2$feature == "Empire\nviews", plot_col, "grey30"), 
             fill = "white", size = 2.5, pch = 21) +
  facet_grid(feature ~ ., scales='free', space='free') +
  labs(y=NULL, x='AMCE', title = "Profile rating") +
  theme_classic() +
  theme(
    strip.placement = 'outside',
    panel.background = element_rect(color='black'),
    legend.position = "bottom",
    legend.title = element_blank(),
    strip.background = element_blank(),
    strip.text = element_text(size = 12),
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(face="bold", size = 10),
    axis.ticks.y = element_blank(),
    axis.text.y = element_blank()) +
  theme(legend.position = "none")


## Combining the choice and the ranking

fullchoice_amce = rbind(amce1, amce2)
ggarrange(conjfig1.1, 
          conjfig1.2,
          widths = c(1, 0.8))
ggsave("conj_comb_amce.png", 
       dpi=500, height=20, width=20, units="cm")



### Conjoint subgroup analysis, MMs----


## Split by imperial nostalgia

plot_col = brewer.pal(5, "BrBG")[5]

dat_cnj$imp_nost_di = ifelse(dat_cnj$imp_att < -0.008258, "Not nostalgic", "Nostalgic") # less than median 
dat_cnj$imp_nost_di = as.factor(dat_cnj$imp_nost_di)

mm3 = cj(MP_chosen ~ MP_empire_views + MP_educ + MP_occup + MP_party + MP_ethnic + MP_gender
          + MP_age + MP_tax_views, 
          data = dat_cnj, id = ~ID, by = ~ imp_nost_di, estimate = "mm", weights = ~wt,
          feature_labels = list("MP_empire_views" = "Empire\nviews",
                                "MP_educ" = "Education", 
                                "MP_occup" = "Occupation",
                                "MP_party" = "Party",
                                "MP_ethnic" = "Ethnicity",
                                "MP_gender" = "Gender", 
                                "MP_age" = "Age",
                                "MP_tax_views" = "Tax\nviews"))

panel1_impnost = mm3 %>% 
  ggplot(aes(level)) +
  geom_vline(xintercept = 0.5, color="grey10") +
  geom_segment(aes(x = lower, xend = upper, y = level, yend = level), 
               color = ifelse(mm3$feature %in% "Empire\nviews", plot_col, "grey30"), 
               size = 2.5, alpha = .8) +
  geom_point(aes(x = estimate, y = level), 
             color = ifelse(mm3$feature == "Empire\nviews", plot_col, "grey30"), 
             fill = "white", size = 2.5, pch = 21) +
  facet_grid(feature ~ imp_nost_di, scales='free', space='free') +
  labs(y=NULL, x='Marginal means') +
  theme_classic() +
  theme(
    strip.placement = 'outside',
    panel.background = element_rect(color='black'),
    legend.position = "bottom",
    legend.title = element_blank(),
    strip.background = element_blank(),
    strip.text = element_text(size = 12),
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(face="bold", size = 10),
    axis.text.y = element_text(face="bold", size = 10)) +
  theme(legend.position = "none")

ggsave("conj_choice_impnost.png", 
       dpi=500, height=20, width=20, units="cm")


## Split by authoritarian values

plot_col = brewer.pal(5, "BrBG")[5]

dat_cnj$auth_di = as.factor(
  ifelse(dat_cnj$auth_val < 0.01393, "Libertarian", "Authoritarian")) # median split 

mm4 = cj(MP_chosen ~ MP_empire_views + MP_educ + MP_occup + MP_party + MP_ethnic + MP_gender
          + MP_age + MP_tax_views, 
          data = dat_cnj, id = ~ID, by = ~ auth_di, estimate = "mm", weights = ~wt,
          feature_labels = list("MP_empire_views" = "Empire\nviews",
                                "MP_educ" = "Education", 
                                "MP_occup" = "Occupation",
                                "MP_party" = "Party",
                                "MP_ethnic" = "Ethnicity",
                                "MP_gender" = "Gender", 
                                "MP_age" = "Age",
                                "MP_tax_views" = "Tax\nviews"))

panel1_auth = mm4 %>% 
  ggplot(aes(level)) +
  geom_vline(xintercept = 0.5, color="grey10") +
  geom_segment(aes(x = lower, xend = upper, y = level, yend = level), 
               color = ifelse(mm4$feature %in% "Empire\nviews", plot_col, "grey30"), 
               size = 2.5, alpha = .8) +
  geom_point(aes(x = estimate, y = level), 
             color = ifelse(mm4$feature == "Empire\nviews", plot_col, "grey30"), 
             fill = "white", size = 2.5, pch = 21) +
  facet_grid(feature ~ auth_di, scales='free', space='free') +
  labs(y=NULL, x='Marginal means') +
  theme_classic() +
  theme(
    strip.placement = 'outside',
    panel.background = element_rect(color='black'),
    legend.position = "bottom",
    legend.title = element_blank(),
    strip.background = element_blank(),
    strip.text = element_text(size = 12),
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(face="bold", size = 10),
    axis.text.y = element_text(face="bold", size = 10)) +
  theme(legend.position = "none")

ggsave("conj_choice_auth.png", 
       dpi=500, height=20, width=20, units="cm")


## Split by national chauvinism

plot_col = brewer.pal(5, "BrBG")[5]

dat_cnj$chauv_di = as.factor(
  ifelse(dat_cnj$nat_chauv < -0.022530, "Cosmopolitan", "Chauvinist")) # median split 

mm5 = cj(MP_chosen ~ MP_empire_views + MP_educ + MP_occup + MP_party + MP_ethnic + MP_gender
          + MP_age + MP_tax_views, 
          data = dat_cnj, id = ~ID, by = ~ chauv_di, estimate = "mm", weights = ~wt,
          feature_labels = list("MP_empire_views" = "Empire\nviews",
                                "MP_educ" = "Education", 
                                "MP_occup" = "Occupation",
                                "MP_party" = "Party",
                                "MP_ethnic" = "Ethnicity",
                                "MP_gender" = "Gender", 
                                "MP_age" = "Age",
                                "MP_tax_views" = "Tax\nviews"))

panel1_chauv = mm5 %>% 
  ggplot(aes(level)) +
  geom_vline(xintercept = 0.5, color="grey10") +
  geom_segment(aes(x = lower, xend = upper, y = level, yend = level), 
               color = ifelse(mm5$feature %in% "Empire\nviews", plot_col, "grey30"), 
               size = 2.5, alpha = .8) +
  geom_point(aes(x = estimate, y = level), 
             color = ifelse(mm5$feature == "Empire\nviews", plot_col, "grey30"), 
             fill = "white", size = 2.5, pch = 21) +
  facet_grid(feature ~ chauv_di, scales='free', space='free') +
  labs(y=NULL, x='Marginal means') +
  theme_classic() +
  theme(
    strip.placement = 'outside',
    panel.background = element_rect(color='black'),
    legend.position = "bottom",
    legend.title = element_blank(),
    strip.background = element_blank(),
    strip.text = element_text(size = 12),
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(face="bold", size = 10),
    axis.text.y = element_text(face="bold", size = 10)) +
  theme(legend.position = "none")

ggsave("conj_choice_chauv.png", 
       dpi=500, height=20, width=20, units="cm")


## split by vote intention - general left v right

plot_col = brewer.pal(5, "BrBG")[5]

dat_cnj = dat_cnj %>% mutate(
  vote_intent_di = case_when(vote_intent == "Cons" ~ "Right party",
                             vote_intent == "Green" ~ "Left party",
                             vote_intent == "Labour" ~ "Left party",
                             vote_intent == "LibDem" ~ "Left party",
                             vote_intent == "Plaid" ~ "Left party",
                             vote_intent == "Reform" ~ "Right party",
                             vote_intent == "SNP" ~ "Left party",
                             vote_intent == "Other" ~ NA,
                             vote_intent == "Not vote" ~ NA))

dat_cnj$vote_intent_di = factor(dat_cnj$vote_intent_di,
                                 levels = c("Right party", "Left party"))

mm8 = cj(MP_chosen ~ MP_empire_views + MP_educ + MP_occup + MP_party + MP_ethnic + MP_gender
          + MP_age + MP_tax_views, 
          data = dat_cnj, id = ~ID, by = ~vote_intent_di, estimate = "mm", weights = ~wt,
          feature_labels = list("MP_empire_views" = "Empire\nviews",
                                "MP_educ" = "Education", 
                                "MP_occup" = "Occupation",
                                "MP_party" = "Party",
                                "MP_ethnic" = "Ethnicity",
                                "MP_gender" = "Gender", 
                                "MP_age" = "Age",
                                "MP_tax_views" = "Tax\nviews"))

mm8 %>% 
  ggplot(aes(level)) +
  geom_vline(xintercept = 0.5, color="grey10") +
  geom_segment(aes(x = lower, xend = upper, y = level, yend = level), 
               color = ifelse(mm8$feature %in% "Empire\nviews", plot_col, "grey30"), 
               size = 2.5, alpha = .8) +
  geom_point(aes(x = estimate, y = level), 
             color = ifelse(mm8$feature == "Empire\nviews", plot_col, "grey30"), 
             fill = "white", size = 2.5, pch = 21) +
  facet_grid(feature~vote_intent_di, scales='free', space='free') +
  labs(y=NULL, x='Marginal means') +
  theme_classic() +
  theme(
    strip.placement = 'outside',
    panel.background = element_rect(color='black'),
    legend.position = "bottom",
    legend.title = element_blank(),
    strip.background= element_blank(),
    strip.text = element_text(size = 12),
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(face="bold", size = 10),
    axis.text.y = element_text(face="bold", size = 10)) +
  theme(legend.position = "none") 

ggsave("conj_choice_vote_int.png", 
       dpi=500, height=20, width=20, units="cm")


## Plots imperial marginal means only for the split samples

mm3 = mm3 %>% 
  filter(feature %in% "Empire\nviews") %>% 
  dplyr::select(-imp_nost_di)
mm3$variable = "Imperial nostalgia"

mm4 = mm4 %>% 
  filter(feature %in% "Empire\nviews") %>% 
  dplyr::select(-auth_di)
mm4$variable = "Authoritarianism"

mm5 = mm5 %>% 
  filter(feature %in% "Empire\nviews") %>% 
  dplyr::select(-chauv_di)
mm5$variable = "National chauvinism"

mm8 = mm8 %>% 
  filter(feature %in% "Empire\nviews") %>% 
  dplyr::select(-vote_intent_di)
mm8$variable = "Voting intention"
mm8$BY = factor(mm8$BY, levels = unique(mm8$BY))

uniform_scale = c(0.3, 0.65) 

nost_split = mm3 %>% 
  ggplot(aes(level)) +
  geom_vline(xintercept = 0.5, color = "grey10") +
  geom_segment(aes(x = lower, xend = upper, y = level, yend = level), 
               size = 2.5, alpha = .8, colour = "#018571") + 
  geom_point(aes(x = estimate, y = level), fill = "white", color = "grey67", 
             size = 2.5, pch = 21) +
  facet_grid(variable ~ BY, scales = 'free', space = 'free') +
  labs(y = NULL, x = NULL, title = NULL) +
  scale_x_continuous(limits = uniform_scale) +
  theme_classic() +
  theme(
    plot.margin = unit(c(.2, .2, .2, .2), "cm"),
    strip.placement = 'outside',
    axis.text.x = element_blank(),
    panel.background = element_rect(color='black'),
    legend.title = element_blank(),
    strip.background = element_blank(),
    strip.text = element_text(size = 12),
    axis.title.x = element_text(face = "bold", size = 10),
    axis.text.y = element_text(face = "bold", size = 10)) +
  theme(legend.position = "none") 

auth_split = mm4 %>% 
  ggplot(aes(level)) +
  geom_vline(xintercept = 0.5, color = "grey10") +
  geom_segment(aes(x = lower, xend = upper, y = level, yend = level), 
               size = 2.5, alpha = .8, colour = "#018571") + 
  geom_point(aes(x = estimate, y = level), fill = "white", color = "grey67", 
             size = 2.5, pch = 21) +
  facet_grid(variable ~ BY, scales = 'free', space = 'free') +
  labs(y = NULL, x = NULL, title = NULL) +
  scale_x_continuous(limits = uniform_scale) +
  theme_classic() +
  theme(
    plot.margin = unit(c(.2, .2, .2, .2), "cm"),
    strip.placement = 'outside',
    axis.text.x = element_blank(),
    panel.background = element_rect(color='black'),
    legend.title = element_blank(),
    strip.background = element_blank(),
    strip.text = element_text(size = 12),
    axis.title.x = element_text(face = "bold", size = 10),
    axis.text.y = element_text(face = "bold", size = 10)) +
  theme(legend.position = "none") 

chauv_split = mm5 %>% 
  ggplot(aes(level)) +
  geom_vline(xintercept = 0.5, color = "grey10") +
  geom_segment(aes(x = lower, xend = upper, y = level, yend = level), 
               size = 2.5, alpha = .8, colour = "#018571") + 
  geom_point(aes(x = estimate, y = level), fill = "white", color = "grey67", 
             size = 2.5, pch = 21) +
  facet_grid(variable ~ BY, scales = 'free', space = 'free') +
  labs(y = NULL, x = NULL, title = NULL) +
  scale_x_continuous(limits = uniform_scale) +
  theme_classic() +
  theme(
    plot.margin = unit(c(.2, .2, .2, .2), "cm"),
    strip.placement = 'outside',
    axis.text.x = element_blank(),
    panel.background = element_rect(color='black'),
    legend.title = element_blank(),
    strip.background = element_blank(),
    strip.text = element_text(size = 12),
    axis.title.x = element_text(face = "bold", size = 10),
    axis.text.y = element_text(face = "bold", size = 10)) +
  theme(legend.position = "none")

vote_split = mm8 %>% 
  ggplot(aes(level)) +
  geom_vline(xintercept = 0.5, color = "grey10") +
  geom_segment(aes(x = lower, xend = upper, y = level, yend = level), 
               size = 2.5, alpha = .8, colour = "#018571") + 
  geom_point(aes(x = estimate, y = level), fill = "white", color = "grey67", 
             size = 2.5, pch = 21) +
  facet_grid(variable ~ BY, scales = 'free', space = 'free') +
  labs(y = NULL, x = NULL, title = NULL) +
  scale_x_continuous(limits = uniform_scale) +
  theme_classic() +
  theme(
    plot.margin = unit(c(.2, .2, .2, .2), "cm"),
    strip.placement = 'outside',
    panel.background = element_rect(color='black'),
    legend.title = element_blank(),
    strip.background = element_blank(),
    strip.text = element_text(size = 12),
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(face = "bold", size = 10),
    axis.text.y = element_text(face = "bold", size = 10)) +
  theme(legend.position = "none")


# Combine the plots and save
combined_plot = ggarrange(
  nost_split, auth_split, chauv_split, vote_split,  
  ncol = 1, nrow = 4, heights = c(1, 1, 1, 1.1),  
  align = "v", vjust = 1, label.x = 0.5  
)
ggsave("mm_subgroup.png", 
       dpi=500, height=18, width=20, units="cm")




###  Conjoint subgroup analysis, AMCEs----


## Split by imperial nostalgia

plot_col = brewer.pal(5, "BrBG")[5]

dat_cnj$imp_nost_di = ifelse(dat_cnj$imp_att < -0.008258, "Not nostalgic", "Nostalgic") # less than median 
dat_cnj$imp_nost_di = as.factor(dat_cnj$imp_nost_di)

amce3 = cj(MP_chosen ~ MP_empire_views + MP_educ + MP_occup + MP_party + MP_ethnic + MP_gender
          + MP_age + MP_tax_views, 
          data = dat_cnj, id = ~ID, by = ~ imp_nost_di, estimate = "amce", weights = ~wt,
          feature_labels = list("MP_empire_views" = "Empire\nviews",
                                "MP_educ" = "Education", 
                                "MP_occup" = "Occupation",
                                "MP_party" = "Party",
                                "MP_ethnic" = "Ethnicity",
                                "MP_gender" = "Gender", 
                                "MP_age" = "Age",
                                "MP_tax_views" = "Tax\nviews"))

panel1_impnost = amce3 %>% 
  ggplot(aes(level)) +
  geom_vline(xintercept = 0, color="grey10") +
  geom_segment(aes(x = lower, xend = upper, y = level, yend = level), 
               color = ifelse(amce3$feature %in% "Empire\nviews", plot_col, "grey30"), 
               size = 2.5, alpha = .8) +
  geom_point(aes(x = estimate, y = level), 
             color = ifelse(amce3$feature == "Empire\nviews", plot_col, "grey30"), 
             fill = "white", size = 2.5, pch = 21) +
  facet_grid(feature ~ imp_nost_di, scales='free', space='free') +
  labs(y=NULL, x='ACME') +
  theme_classic() +
  theme(
    strip.placement = 'outside',
    panel.background = element_rect(color='black'),
    legend.position = "bottom",
    legend.title = element_blank(),
    strip.background = element_blank(),
    strip.text = element_text(size = 12),
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(face="bold", size = 10),
    axis.text.y = element_text(face="bold", size = 10)) +
  theme(legend.position = "none")


## Split by authoritarian values

plot_col = brewer.pal(5, "BrBG")[5]

dat_cnj$auth_di = as.factor(
  ifelse(dat_cnj$auth_val < 0.01393, "Libertarian", "Authoritarian")) # median split 

amce4 = cj(MP_chosen ~ MP_empire_views + MP_educ + MP_occup + MP_party + MP_ethnic + MP_gender
          + MP_age + MP_tax_views, 
          data = dat_cnj, id = ~ID, by = ~ auth_di, estimate = "amce", weights = ~wt,
          feature_labels = list("MP_empire_views" = "Empire\nviews",
                                "MP_educ" = "Education", 
                                "MP_occup" = "Occupation",
                                "MP_party" = "Party",
                                "MP_ethnic" = "Ethnicity",
                                "MP_gender" = "Gender", 
                                "MP_age" = "Age",
                                "MP_tax_views" = "Tax\nviews"))

panel1_auth = amce4 %>% 
  ggplot(aes(level)) +
  geom_vline(xintercept = 0, color="grey10") +
  geom_segment(aes(x = lower, xend = upper, y = level, yend = level), 
               color = ifelse(amce4$feature %in% "Empire\nviews", plot_col, "grey30"), 
               size = 2.5, alpha = .8) +
  geom_point(aes(x = estimate, y = level), 
             color = ifelse(amce4$feature == "Empire\nviews", plot_col, "grey30"), 
             fill = "white", size = 2.5, pch = 21) +
  facet_grid(feature ~ auth_di, scales='free', space='free') +
  labs(y=NULL, x='ACME') +
  theme_classic() +
  theme(
    strip.placement = 'outside',
    panel.background = element_rect(color='black'),
    legend.position = "bottom",
    legend.title = element_blank(),
    strip.background = element_blank(),
    strip.text = element_text(size = 12),
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(face="bold", size = 10),
    axis.text.y = element_text(face="bold", size = 10)) +
  theme(legend.position = "none")


## Split by national chauvinism

plot_col = brewer.pal(5, "BrBG")[5]

dat_cnj$chauv_di = as.factor(
  ifelse(dat_cnj$nat_chauv < -0.022530, "Cosmopolitan", "Chauvinist")) # median split 

amce5 = cj(MP_chosen ~ MP_empire_views + MP_educ + MP_occup + MP_party + MP_ethnic + MP_gender
          + MP_age + MP_tax_views, 
          data = dat_cnj, id = ~ID, by = ~ chauv_di, estimate = "amce", weights = ~wt,
          feature_labels = list("MP_empire_views" = "Empire\nviews",
                                "MP_educ" = "Education", 
                                "MP_occup" = "Occupation",
                                "MP_party" = "Party",
                                "MP_ethnic" = "Ethnicity",
                                "MP_gender" = "Gender", 
                                "MP_age" = "Age",
                                "MP_tax_views" = "Tax\nviews"))

panel1_chauv = amce5 %>% 
  ggplot(aes(level)) +
  geom_vline(xintercept = 0, color="grey10") +
  geom_segment(aes(x = lower, xend = upper, y = level, yend = level), 
               color = ifelse(amce5$feature %in% "Empire\nviews", plot_col, "grey30"), 
               size = 2.5, alpha = .8) +
  geom_point(aes(x = estimate, y = level), 
             color = ifelse(amce5$feature == "Empire\nviews", plot_col, "grey30"), 
             fill = "white", size = 2.5, pch = 21) +
  facet_grid(feature ~ chauv_di, scales='free', space='free') +
  labs(y=NULL, x='ACME') +
  theme_classic() +
  theme(
    strip.placement = 'outside',
    panel.background = element_rect(color='black'),
    legend.position = "bottom",
    legend.title = element_blank(),
    strip.background = element_blank(),
    strip.text = element_text(size = 12),
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(face="bold", size = 10),
    axis.text.y = element_text(face="bold", size = 10)) +
  theme(legend.position = "none")


## split by vote intention - general left v right

plot_col = brewer.pal(5, "BrBG")[5]

dat_cnj = dat_cnj %>% mutate(
  vote_intent_di = case_when(vote_intent == "Cons" ~ "Right party",
                             vote_intent == "Green" ~ "Left party",
                             vote_intent == "Labour" ~ "Left party",
                             vote_intent == "LibDem" ~ "Left party",
                             vote_intent == "Plaid" ~ "Left party",
                             vote_intent == "Reform" ~ "Right party",
                             vote_intent == "SNP" ~ "Left party",
                             vote_intent == "Other" ~ NA,
                             vote_intent == "Not vote" ~ NA))

dat_cnj$vote_intent_di = factor(dat_cnj$vote_intent_di,
                                 levels = c("Right party", "Left party"))

amce8 = cj(MP_chosen ~ MP_empire_views + MP_educ + MP_occup + MP_party + MP_ethnic + MP_gender
          + MP_age + MP_tax_views, 
          data = dat_cnj, id = ~ID, by = ~vote_intent_di, estimate = "amce", weights = ~wt,
          feature_labels = list("MP_empire_views" = "Empire\nviews",
                                "MP_educ" = "Education", 
                                "MP_occup" = "Occupation",
                                "MP_party" = "Party",
                                "MP_ethnic" = "Ethnicity",
                                "MP_gender" = "Gender", 
                                "MP_age" = "Age",
                                "MP_tax_views" = "Tax\nviews"))

amce8 %>% 
  ggplot(aes(level)) +
  geom_vline(xintercept = 0, color="grey10") +
  geom_segment(aes(x = lower, xend = upper, y = level, yend = level), 
               color = ifelse(amce8$feature %in% "Empire\nviews", plot_col, "grey30"), 
               size = 2.5, alpha = .8) +
  geom_point(aes(x = estimate, y = level), 
             color = ifelse(amce8$feature == "Empire\nviews", plot_col, "grey30"), 
             fill = "white", size = 2.5, pch = 21) +
  facet_grid(feature~vote_intent_di, scales='free', space='free') +
  labs(y=NULL, x='ACME') +
  theme_classic() +
  theme(
    strip.placement = 'outside',
    panel.background = element_rect(color='black'),
    legend.position = "bottom",
    legend.title = element_blank(),
    strip.background= element_blank(),
    strip.text = element_text(size = 12),
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(face="bold", size = 10),
    axis.text.y = element_text(face="bold", size = 10)) +
  theme(legend.position = "none") 


## Plots imperial marginal means only for the split samples

amce3 = amce3 %>% 
  filter(feature %in% "Empire\nviews") %>% 
  dplyr::select(-imp_nost_di)
amce3$variable = "Imperial nostalgia"

amce4 = amce4 %>% 
  filter(feature %in% "Empire\nviews") %>% 
  dplyr::select(-auth_di)
amce4$variable = "Authoritarianism"

amce5 = amce5 %>% 
  filter(feature %in% "Empire\nviews") %>% 
  dplyr::select(-chauv_di)
amce5$variable = "National chauvinism"

amce8 = amce8 %>% 
  filter(feature %in% "Empire\nviews") %>% 
  dplyr::select(-vote_intent_di)
amce8$variable = "Voting intention"
amce8$BY = factor(amce8$BY, levels = unique(amce8$BY))

uniform_scale = c(-0.25, 0.25) 

nost_split = amce3 %>% 
  ggplot(aes(level)) +
  geom_vline(xintercept = 0, color = "grey10") +
  geom_segment(aes(x = lower, xend = upper, y = level, yend = level), 
               size = 2.5, alpha = .8, colour = "#018571") + 
  geom_point(aes(x = estimate, y = level), fill = "white", color = "grey67", 
             size = 2.5, pch = 21) +
  facet_grid(variable ~ BY, scales = 'free', space = 'free') +
  labs(y = NULL, x = NULL, title = NULL) +
  scale_x_continuous(limits = uniform_scale) +
  theme_classic() +
  theme(
    plot.margin = unit(c(.2, .2, .2, .2), "cm"),
    strip.placement = 'outside',
    axis.text.x = element_blank(),
    panel.background = element_rect(color='black'),
    legend.title = element_blank(),
    strip.background = element_blank(),
    strip.text = element_text(size = 12),
    axis.title.x = element_text(face = "bold", size = 10),
    axis.text.y = element_text(face = "bold", size = 10)) +
  theme(legend.position = "none") 

auth_split = amce4 %>% 
  ggplot(aes(level)) +
  geom_vline(xintercept = 0, color = "grey10") +
  geom_segment(aes(x = lower, xend = upper, y = level, yend = level), 
               size = 2.5, alpha = .8, colour = "#018571") + 
  geom_point(aes(x = estimate, y = level), fill = "white", color = "grey67", 
             size = 2.5, pch = 21) +
  facet_grid(variable ~ BY, scales = 'free', space = 'free') +
  labs(y = NULL, x = NULL, title = NULL) +
  scale_x_continuous(limits = uniform_scale) +
  theme_classic() +
  theme(
    plot.margin = unit(c(.2, .2, .2, .2), "cm"),
    strip.placement = 'outside',
    axis.text.x = element_blank(),
    panel.background = element_rect(color='black'),
    legend.title = element_blank(),
    strip.background = element_blank(),
    strip.text = element_text(size = 12),
    axis.title.x = element_text(face = "bold", size = 10),
    axis.text.y = element_text(face = "bold", size = 10)) +
  theme(legend.position = "none") 

chauv_split = amce5 %>% 
  ggplot(aes(level)) +
  geom_vline(xintercept = 0, color = "grey10") +
  geom_segment(aes(x = lower, xend = upper, y = level, yend = level), 
               size = 2.5, alpha = .8, colour = "#018571") + 
  geom_point(aes(x = estimate, y = level), fill = "white", color = "grey67", 
             size = 2.5, pch = 21) +
  facet_grid(variable ~ BY, scales = 'free', space = 'free') +
  labs(y = NULL, x = NULL, title = NULL) +
  scale_x_continuous(limits = uniform_scale) +
  theme_classic() +
  theme(
    plot.margin = unit(c(.2, .2, .2, .2), "cm"),
    strip.placement = 'outside',
    axis.text.x = element_blank(),
    panel.background = element_rect(color='black'),
    legend.title = element_blank(),
    strip.background = element_blank(),
    strip.text = element_text(size = 12),
    axis.title.x = element_text(face = "bold", size = 10),
    axis.text.y = element_text(face = "bold", size = 10)) +
  theme(legend.position = "none")

vote_split = amce8 %>% 
  ggplot(aes(level)) +
  geom_vline(xintercept = 0, color = "grey10") +
  geom_segment(aes(x = lower, xend = upper, y = level, yend = level), 
               size = 2.5, alpha = .8, colour = "#018571") + 
  geom_point(aes(x = estimate, y = level), fill = "white", color = "grey67", 
             size = 2.5, pch = 21) +
  facet_grid(variable ~ BY, scales = 'free', space = 'free') +
  labs(y = NULL, x = NULL, title = NULL) +
  scale_x_continuous(limits = uniform_scale) +
  theme_classic() +
  theme(
    plot.margin = unit(c(.2, .2, .2, .2), "cm"),
    strip.placement = 'outside',
    panel.background = element_rect(color='black'),
    legend.title = element_blank(),
    strip.background = element_blank(),
    strip.text = element_text(size = 12),
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(face = "bold", size = 10),
    axis.text.y = element_text(face = "bold", size = 10)) +
  theme(legend.position = "none")


# Combine the plots and save
combined_plot = ggarrange(
  nost_split, auth_split, chauv_split, vote_split,  
  ncol = 1, nrow = 4, heights = c(1, 1, 1, 1.1),  
  align = "v", vjust = 1, label.x = 0.5  
)
ggsave("amce_subgroup.png", 
       dpi=500, height=18, width=20, units="cm")



### Difference in marginal means----

# imperial nostalgia

dat_cnj$imp_nost_di = as.factor(
  ifelse(dat_cnj$imp_att < -0.008258, "Not nostalgic", "Nostalgic")) # less than median 

mm_imp_nost = cj(MP_chosen ~ MP_empire_views + MP_educ + MP_occup + MP_party 
                  + MP_ethnic + MP_gender + MP_age + MP_tax_views, 
                  data = dat_cnj, id = ~ID, by = ~ imp_nost_di, estimate = "mm_differences", 
                  weights = ~wt, 
                  feature_labels = list("MP_empire_views" = "Empire\nviews",
                                        "MP_educ" = "Education", 
                                        "MP_occup" = "Occupation", 
                                        "MP_party" = "Party",
                                        "MP_ethnic" = "Ethnicity",
                                        "MP_gender" = "Gender", 
                                        "MP_age" = "Age",
                                        "MP_tax_views" = "Tax\nviews"))

mm_imp_nost = mm_imp_nost %>% 
  filter(feature %in% "Empire\nviews") %>% 
  dplyr::select(-imp_nost_di)
mm_imp_nost$variable = "Imperial nostalgia"

# authoritarianism

dat_cnj$auth_di = as.factor(
  ifelse(dat_cnj$auth_val < 0.01393, "Libertarian", "Authoritarian")) # median split 

mm_auth = cj(MP_chosen ~ MP_empire_views + MP_educ + MP_occup + MP_party 
              + MP_ethnic + MP_gender + MP_age + MP_tax_views, 
              data = dat_cnj, id = ~ID, by = ~ auth_di, estimate = "mm_differences", 
              weights = ~wt,
              feature_labels = list("MP_empire_views" = "Empire\nviews",
                                    "MP_educ" = "Education", 
                                    "MP_occup" = "Occupation",
                                    "MP_party" = "Party",
                                    "MP_ethnic" = "Ethnicity",
                                    "MP_gender" = "Gender", 
                                    "MP_age" = "Age",
                                    "MP_tax_views" = "Tax\nviews"))

mm_auth = mm_auth %>% 
  filter(feature %in% "Empire\nviews") %>% 
  dplyr::select(-auth_di)
mm_auth$variable = "Authoritarianism"

# national chauvinism

dat_cnj$chauv_di = as.factor(
  ifelse(dat_cnj$nat_chauv < -0.022530, "Cosmopolitan", "Chauvinist")) # median split 

mm_chauv = cj(MP_chosen ~ MP_empire_views + MP_educ + MP_occup + MP_party 
               + MP_ethnic + MP_gender + MP_age + MP_tax_views, 
               data = dat_cnj, id = ~ID, by = ~ chauv_di, estimate = "mm_differences", 
               weights = ~wt,
               feature_labels = list("MP_empire_views" = "Empire\nviews",
                                     "MP_educ" = "Education", 
                                     "MP_occup" = "Occupation",
                                     "MP_party" = "Party",
                                     "MP_ethnic" = "Ethnicity",
                                     "MP_gender" = "Gender", 
                                     "MP_age" = "Age",
                                     "MP_tax_views" = "Tax\nviews"))

mm_chauv = mm_chauv %>% 
  filter(feature %in% "Empire\nviews") %>% 
  dplyr::select(-chauv_di)
mm_chauv$variable = "National chauvinism"

# vote intention - left v right

dat_cnj = dat_cnj %>% mutate(
  vote_intent_di = case_when(vote_intent == "Cons" ~ "Right party",
                             vote_intent == "Green" ~ "Left party",
                             vote_intent == "Labour" ~ "Left party",
                             vote_intent == "LibDem" ~ "Left party",
                             vote_intent == "Plaid" ~ "Left party",
                             vote_intent == "Reform" ~ "Right party",
                             vote_intent == "SNP" ~ "Left party",
                             vote_intent == "Other" ~ NA,
                             vote_intent == "Not vote" ~ NA))

dat_cnj$vote_intent_di = factor(dat_cnj$vote_intent_di,
                                 levels = c("Right party", "Left party"))

mm_vote = cj(MP_chosen ~ MP_empire_views + MP_educ + MP_occup + MP_party 
              + MP_ethnic + MP_gender + MP_age + MP_tax_views, 
              data = dat_cnj, id = ~ID, by = ~vote_intent_di, estimate = "mm_differences", 
              weights = ~wt, 
              feature_labels = list("MP_empire_views" = "Empire\nviews",
                                    "MP_educ" = "Education", 
                                    "MP_occup" = "Occupation",
                                    "MP_party" = "Party",
                                    "MP_ethnic" = "Ethnicity",
                                    "MP_gender" = "Gender", 
                                    "MP_age" = "Age",
                                    "MP_tax_views" = "Tax\nviews"))

mm_vote = mm_vote %>% 
  filter(feature %in% "Empire\nviews") %>% 
  dplyr::select(-vote_intent_di)
mm_vote$variable = "Voting intention"

# plot

mm_all_covar = bind_rows(mm_imp_nost, mm_auth, mm_chauv, mm_vote)

mm_all_covar %>% ggplot(aes(level)) +
  geom_vline(xintercept = 0, color="grey10") +
  geom_segment(aes(x=lower, xend=upper, y=level, yend=level), size=3, alpha=.8, 
               colour = "#018571") + 
  geom_point(aes(x=estimate, y=level), fill="white", color="grey67", size=2.5,  pch=21)+
  facet_grid(variable ~ ., scales='free', space='free') +
  #coord_flip() +
  labs(y=NULL, x='', title = "Differences in marginal means") +
  theme_classic() +
  theme(
    strip.placement = 'outside',
    panel.background = element_rect(color='black'),
    legend.position = "bottom",
    legend.title = element_blank(),
    strip.background = element_blank(),
    strip.text = element_text(size = 12),
    axis.text.x = element_text(size = 10),
    axis.title.x = element_text(face="bold", size = 10),
    axis.text.y = element_text(face="bold", size = 10)) +
  theme(legend.position = "none") +
  geom_text(
    aes(label = sprintf("%0.2f", round(estimate, digits = 3)), y = level, x = estimate),
    colour = "#018571",
    position = position_dodge(width = 0.9), size = 4, vjust = 1.5, fontface = "bold")

ggsave("mmdiffs_subgroup.png", 
       dpi=500, height=18, width=20, units="cm")



## Combined plot for marginal means by subgroup and differences in MMs


# Fix variable order to match
mm_all_covar$variable = factor(mm_all_covar$variable, 
                                levels = c("Imperial nostalgia", "Authoritarianism", 
                                           "National chauvinism", "Voting intention"))

# x lims
subgrp_scale = c(0.3, 0.65) 
diff_scale = c(-0.2, 0.25) 

# Create individual plots for each variable and subgroup
# Row 1: Imperial nostalgia
nost_left = mm3 %>% 
  filter(BY == "Not nostalgic") %>%
  ggplot(aes(level)) +
  geom_vline(xintercept = 0.5, color = "grey10") +
  geom_segment(aes(x = lower, xend = upper, y = level, yend = level), 
               size = 2.5, alpha = .8, colour = "#018571") + 
  geom_point(aes(x = estimate, y = level), fill = "white", color = "grey67", 
             size = 2.5, pch = 21) +
  labs(y = "Imperial nostalgia", x = NULL, title = "Not nostalgic") +
  scale_x_continuous(limits = subgrp_scale) +
  theme_classic() +
  theme(
    plot.margin = unit(c(.2, .1, .2, .2), "cm"),
    panel.background = element_rect(color='black'),
    plot.title = element_text(size = 12, hjust = 0.5),
    axis.text.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.y = element_text(face = "bold", size = 10))

nost_right = mm3 %>% 
  filter(BY == "Nostalgic") %>%
  ggplot(aes(level)) +
  geom_vline(xintercept = 0.5, color = "grey10") +
  geom_segment(aes(x = lower, xend = upper, y = level, yend = level), 
               size = 2.5, alpha = .8, colour = "#018571") + 
  geom_point(aes(x = estimate, y = level), fill = "white", color = "grey67", 
             size = 2.5, pch = 21) +
  labs(y = NULL, x = NULL, title = "Nostalgic") +
  scale_x_continuous(limits = subgrp_scale) +
  theme_classic() +
  theme(
    plot.margin = unit(c(.2, .1, .2, .1), "cm"),
    panel.background = element_rect(color='black'),
    plot.title = element_text(size = 12, hjust = 0.5),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank())

nost_diff = mm_all_covar %>% 
  filter(variable == "Imperial nostalgia") %>%
  ggplot(aes(level)) +
  geom_vline(xintercept = 0, color="grey10") +
  geom_segment(aes(x = lower, xend = upper, y = level, yend = level), 
               size = 2.5, alpha = .8, colour = "#018571") + 
  geom_point(aes(x=estimate, y=level), fill="white", color="grey67", size=2.5, pch=21)+
  labs(y=NULL, x='', title = "Difference (left - right)") +
  scale_x_continuous(limits = diff_scale) +
  theme_classic() +
  theme(
    plot.margin = unit(c(.2, .1, -.23, .1), "cm"),
    panel.background = element_rect(color='black'),
    plot.title = element_text(size = 12, hjust = 0.5),
    axis.title.y = element_text(size = 12, angle = 90),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank())

# Row 2: Authoritarianism
auth_left = mm4 %>% 
  filter(BY == "Libertarian") %>%
  ggplot(aes(level)) +
  geom_vline(xintercept = 0.5, color = "grey10") +
  geom_segment(aes(x = lower, xend = upper, y = level, yend = level), 
               size = 2.5, alpha = .8, colour = "#018571") + 
  geom_point(aes(x = estimate, y = level), fill = "white", color = "grey67", 
             size = 2.5, pch = 21) +
  labs(y = "Authoritarianism", x = NULL, title = "Libertarian") +
  scale_x_continuous(limits = subgrp_scale) +
  theme_classic() +
  theme(
    plot.margin = unit(c(.2, .1, .2, .2), "cm"),
    panel.background = element_rect(color='black'),
    plot.title = element_text(size = 12, hjust = 0.5),
    axis.text.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.y = element_text(face = "bold", size = 10))

auth_right = mm4 %>% 
  filter(BY == "Authoritarian") %>%
  ggplot(aes(level)) +
  geom_vline(xintercept = 0.5, color = "grey10") +
  geom_segment(aes(x = lower, xend = upper, y = level, yend = level), 
               size = 2.5, alpha = .8, colour = "#018571") + 
  geom_point(aes(x = estimate, y = level), fill = "white", color = "grey67", 
             size = 2.5, pch = 21) +
  labs(y = NULL, x = NULL, title = "Authoritarian") +
  scale_x_continuous(limits = subgrp_scale) +
  theme_classic() +
  theme(
    plot.margin = unit(c(.2, .1, .2, .1), "cm"),
    panel.background = element_rect(color='black'),
    plot.title = element_text(size = 12, hjust = 0.5),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank())

auth_diff = mm_all_covar %>% 
  filter(variable == "Authoritarianism") %>%
  ggplot(aes(level)) +
  geom_vline(xintercept = 0, color="grey10") +
  geom_segment(aes(x=lower, xend=upper, y=level, yend=level), size=2.5, alpha=.8, 
               colour = "#018571") + 
  geom_point(aes(x=estimate, y=level), fill="white", color="grey67", size=2.5, pch=21)+
  labs(y=NULL, x='', title = "Difference (left - right)") +
  scale_x_continuous(limits = diff_scale) +
  theme_classic() +
  theme(
    plot.margin = unit(c(.2, .1, -.23, .1), "cm"),
    panel.background = element_rect(color='black'),
    plot.title = element_text(size = 12, hjust = 0.5),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.title.y = element_text(size = 12, angle = 90),
    axis.ticks.y = element_blank()) 

# Row 3: National chauvinism
chauv_left = mm5 %>% 
  filter(BY == "Cosmopolitan") %>%
  ggplot(aes(level)) +
  geom_vline(xintercept = 0.5, color = "grey10") +
  geom_segment(aes(x = lower, xend = upper, y = level, yend = level), 
               size = 2.5, alpha = .8, colour = "#018571") + 
  geom_point(aes(x = estimate, y = level), fill = "white", color = "grey67", 
             size = 2.5, pch = 21) +
  labs(y = "National chauvinism", x = NULL, title = "Cosmopolitan") +
  scale_x_continuous(limits = subgrp_scale) +
  theme_classic() +
  theme(
    plot.margin = unit(c(.2, .1, .2, .1), "cm"),
    panel.background = element_rect(color='black'),
    plot.title = element_text(size = 12, hjust = 0.5),
    axis.text.x = element_blank(),
    axis.title.y = element_blank(),
    axis.text.y = element_text(face = "bold", size = 10))

chauv_right = mm5 %>% 
  filter(BY == "Chauvinist") %>%
  ggplot(aes(level)) +
  geom_vline(xintercept = 0.5, color = "grey10") +
  geom_segment(aes(x = lower, xend = upper, y = level, yend = level), 
               size = 2.5, alpha = .8, colour = "#018571") + 
  geom_point(aes(x = estimate, y = level), fill = "white", color = "grey67", 
             size = 2.5, pch = 21) +
  labs(y = NULL, x = NULL, title = "Chauvinist") +
  scale_x_continuous(limits = subgrp_scale) +
  theme_classic() +
  theme(
    plot.margin = unit(c(.2, .1, .2, .1), "cm"),
    panel.background = element_rect(color='black'),
    plot.title = element_text(size = 12, hjust = 0.5),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank())

chauv_diff = mm_all_covar %>% 
  filter(variable == "National chauvinism") %>%
  ggplot(aes(level)) +
  geom_vline(xintercept = 0, color="grey10") +
  geom_segment(aes(x=lower, xend=upper, y=level, yend=level), size=2.5, alpha=.8, 
               colour = "#018571") + 
  geom_point(aes(x=estimate, y=level), fill="white", color="grey67", size=2.5, pch=21)+
  labs(y=NULL, x='', title = "Difference (left - right)") +
  scale_x_continuous(limits = diff_scale) +
  theme_classic() +
  theme(
    plot.margin = unit(c(.2, .1, -.23, .1), "cm"),
    panel.background = element_rect(color='black'),
    plot.title = element_text(size = 12, hjust = 0.5),
    axis.text.x = element_blank(),
    axis.text.y = element_blank(),
    axis.title.y = element_text(size = 12, angle = 90),
    axis.ticks.y = element_blank()) 

# Row 4: Voting intention
vote_left = mm8 %>% 
  filter(BY == "Left party") %>%
  ggplot(aes(level)) +
  geom_vline(xintercept = 0.5, color = "grey10") +
  geom_segment(aes(x = lower, xend = upper, y = level, yend = level), 
               size = 2.5, alpha = .8, colour = "#018571") + 
  geom_point(aes(x = estimate, y = level), fill = "white", color = "grey67", 
             size = 2.5, pch = 21) +
  labs(y = "Voting intention", x = NULL, title = "Left party") +
  scale_x_continuous(limits = subgrp_scale) +
  theme_classic() +
  theme(
    plot.margin = unit(c(.2, .1, .2, .2), "cm"),
    panel.background = element_rect(color='black'),
    plot.title = element_text(size = 12, hjust = 0.5),
    axis.text.x = element_text(size = 10),
    axis.title.y = element_blank(),
    axis.text.y = element_text(face = "bold", size = 10))

vote_right = mm8 %>% 
  filter(BY == "Right party") %>%
  ggplot(aes(level)) +
  geom_vline(xintercept = 0.5, color = "grey10") +
  geom_segment(aes(x = lower, xend = upper, y = level, yend = level), 
               size = 2.5, alpha = .8, colour = "#018571") + 
  geom_point(aes(x = estimate, y = level), fill = "white", color = "grey67", 
             size = 2.5, pch = 21) +
  labs(y = NULL, x = NULL, title = "Right party") +
  scale_x_continuous(limits = subgrp_scale) +
  theme_classic() +
  theme(
    plot.margin = unit(c(.2, .1, .2, .1), "cm"),
    panel.background = element_rect(color='black'),
    plot.title = element_text(size = 12, hjust = 0.5),
    axis.text.x = element_text(size = 10),
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank())

vote_diff = mm_all_covar %>% 
  filter(variable == "Voting intention") %>%
  ggplot(aes(level)) +
  geom_vline(xintercept = 0, color="grey10") +
  geom_segment(aes(x=lower, xend=upper, y=level, yend=level), size=2.5, alpha=.8, 
               colour = "#018571") + 
  geom_point(aes(x=estimate, y=level), fill="white", color="grey67", size=2.5, pch=21)+
  labs(y=NULL, x='', title = "Difference (left - right)") +
  scale_x_continuous(limits = diff_scale) +
  theme_classic() +
  theme(
    plot.margin = unit(c(.2, .1, -.23, .1), "cm"),
    panel.background = element_rect(color='black'),
    plot.title = element_text(size = 12, hjust = 0.5),
    axis.text.x = element_text(size = 10),
    axis.text.y = element_blank(),
    axis.title.y = element_text(size = 12, angle = 90),
    axis.ticks.y = element_blank()) 

# Arrange in 4×3 grid 
combined_final_plot = ggarrange(
  nost_left, nost_right, nost_diff,
  auth_left, auth_right, auth_diff,
  chauv_left, chauv_right, chauv_diff,
  vote_left, vote_right, vote_diff,
  ncol = 3, nrow = 4,
  widths = c(1.4, 1, 1),
  heights = c(1, 1, 1, 1)
)

# Save the combined plot
ggsave("combined_subgroup_analysis.png", 
       plot = combined_final_plot,
       dpi = 500, height = 18, width = 24, units = "cm")



### Direct analysis of conjoint----


# extract pro- vs anti-empire choices and reshape

cand_choic_emp_di = dat_cnj %>%
  group_by(ID, replicate) %>%
  filter(n() == 2,  # Ensure we have complete pairs
         sum(MP_empire_views == "Civilising effect") == 1 & 
           sum(MP_empire_views == "Atrocities") == 1) %>%
  summarise(
    wt = first(wt),
    imp_att = first(imp_att),
    auth_val = first(auth_val),
    nat_chauv = first(nat_chauv),
    vote_intent_di = first(vote_intent_di),
    sup_Cons = first(sup_Cons),
    sup_Lab = first(sup_Lab),
    sup_Reform = first(sup_Reform),
    sup_LibDem = first(sup_LibDem),
    sup_Green = first(sup_Green),
    pro_empire_chosen = sum((MP_empire_views == "Civilising effect") * MP_chosen),
    .groups = 'drop'
  )

# fit lms with clustered SEs

ch_mod_1 = lm(pro_empire_chosen ~ imp_att + auth_val + sup_Cons + sup_Lab 
              + sup_LibDem + sup_Reform + sup_Green, 
              data = cand_choic_emp_di, weights = wt)
summary(ch_mod_1)
ch_mod_1_cl = coeftest(ch_mod_1, vcov = vcovCL, cluster = ~ID)
ch_mod_1_cl

ch_mod_2 = lm(pro_empire_chosen ~ imp_att + auth_val + vote_intent_di, 
              data = cand_choic_emp_di, weights = wt)
summary(ch_mod_2)
ch_mod_2_cl = coeftest(ch_mod_2, vcov = vcovCL, cluster = ~ID)
ch_mod_2_cl

# table
texreg(list(ch_mod_1, ch_mod_2),
       file = "cand_choice_model.tex",
       override.se = list(ch_mod_1_cl[, "Std. Error"], ch_mod_2_cl[, "Std. Error"]),
       override.pvalues = list(ch_mod_1_cl[, "Pr(>|t|)"], ch_mod_2_cl[, "Pr(>|t|)"]),
       custom.coef.names = c("Intercept", 
                             "Imperial nostalgia", 
                             "Authoritarian values",
                             "Conservative party evaluations", 
                             "Labour party evaluations", 
                             "Liberal Democrat evaluations", 
                             "Reform party evaluations", 
                             "Green party evaluations",
                             "Vote intention for left (vs. right) party"),
       caption = "Conjoint Analysis: Choice of Pro-Empire vs Anti-Empire Candidate",
       label = "tab:empire_choice",
       custom.note = "Clustered standard errors (by respondent ID) in parentheses. %stars.",
       stars = c(0.05), dcolumn = TRUE, booktabs = TRUE, leading.zero = FALSE, digits = 2)


